| title: “Week 3: Use R as GIS” |
| output: github_document |
knitr::opts_chunk$set(fig.width=12, fig.height=8, fig.path='Figs/', warning=FALSE, message=FALSE)
## [[1]]
## [1] "sp" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[2]]
## [1] "rgdal" "sp" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"
##
## [[3]]
## [1] "rgeos" "rgdal" "sp" "stats" "graphics"
## [6] "grDevices" "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "maptools" "rgeos" "rgdal" "sp" "stats"
## [6] "graphics" "grDevices" "utils" "datasets" "methods"
## [11] "base"
##
## [[5]]
## [1] "classInt" "maptools" "rgeos" "rgdal" "sp"
## [6] "stats" "graphics" "grDevices" "utils" "datasets"
## [11] "methods" "base"
##
## [[6]]
## [1] "RColorBrewer" "classInt" "maptools" "rgeos"
## [5] "rgdal" "sp" "stats" "graphics"
## [9] "grDevices" "utils" "datasets" "methods"
## [13] "base"
##
## [[7]]
## [1] "GISTools" "MASS" "RColorBrewer" "classInt"
## [5] "maptools" "rgeos" "rgdal" "sp"
## [9] "stats" "graphics" "grDevices" "utils"
## [13] "datasets" "methods" "base"
##
## [[8]]
## [1] "maps" "GISTools" "MASS" "RColorBrewer"
## [5] "classInt" "maptools" "rgeos" "rgdal"
## [9] "sp" "stats" "graphics" "grDevices"
## [13] "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "raster" "maps" "GISTools" "MASS"
## [5] "RColorBrewer" "classInt" "maptools" "rgeos"
## [9] "rgdal" "sp" "stats" "graphics"
## [13] "grDevices" "utils" "datasets" "methods"
## [17] "base"
##
## [[10]]
## [1] "ggmap" "ggplot2" "raster" "maps"
## [5] "GISTools" "MASS" "RColorBrewer" "classInt"
## [9] "maptools" "rgeos" "rgdal" "sp"
## [13] "stats" "graphics" "grDevices" "utils"
## [17] "datasets" "methods" "base"
| Without attributes | With attributes | |
|---|---|---|
| Points | SpatialPoints | SpatialPointsDataFrame |
| Lines | SpatialLines | SpatialLinesDataFrame |
| Polygons | SpatialPolygons | SpatialPolygonsDataFrame |
| Raster | SpatialGrid | SpatialGridDataFrame |
| Raster | SpatialPixels | SpatialPixelsDataFrame |
LubbockBlock<-readShapePoly("Data/LubbockBlockNew.shp") #read polygon shapefile
class(LubbockBlock)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
HouseLocation<-read.csv("Data/HouseLatLon.csv") #read GPS data
class(HouseLocation)
## [1] "data.frame"
coordinates(HouseLocation)<-c('Lon', 'Lat')
class(HouseLocation)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
cropland<-raster("Data/Lubbock_CDL_2013_USDA.tif")
class(cropland)
## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"
tmin <- getData("worldclim", var = "tmin", res = 10) # this will download
class(tmin)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"
LubbockBlock<-readOGR("./Data", "LubbockBlockNew") #read polygon shapefile
## OGR data source with driver: ESRI Shapefile
## Source: "./Data", layer: "LubbockBlockNew"
## with 167 features
## It has 69 fields
class(LubbockBlock)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
LubbockBlock<-readShapePoly("Data/LubbockBlockNew.shp") #read polygon shapefile
plot(LubbockBlock,axes=TRUE, col=alpha("gray70", 0.6)) #plot Lubbock block shapefile
#add title, scalebar, north arrow, and legend
HouseLocation<-read.csv("Data/HouseLatLon.csv") #read GPS data
price<-HouseLocation$TotalPrice
nclr<-5
priceclr<-brewer.pal(nclr, "Reds")
class<-classIntervals(price, nclr, style="quantile")
clocode<-findColours(class, priceclr)
points(HouseLocation$Lon, HouseLocation$Lat, pch=19, col=clocode, cex=0.5) #add houses on top of Lubbock block shapefile
title(main="Houses on Sale in Lubbock, 2014")
legend(-101.95, 33.65, legend=names(attr(clocode, "table")), fill =attr(clocode, "palette"), cex=0.5, bty="n")
#map.scale(x=-101.85, y=33.49,0.001,"Miles",4,0.5,sfcol='red')
north.arrow(xb=-101.95, yb=33.65, len=0.005, lab="N")
#plot raster
plot(cropland)
#plot raster stack
tmin <- getData("worldclim", var = "tmin", res = 10) # this will download
plot(tmin)
#project a vector
boudary=readShapePoly('Data/boundary');
proj4string(boudary) <-CRS("+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj4string(boudary)
## [1] "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
boudaryProj<-spTransform(boudary, CRS("+init=epsg:3857"))
proj4string(boudaryProj)
## [1] "+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
#project a raster
proj4string(cropland)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(cropland)
aea <- CRS("+init=ESRI:102003") #Albert equal area
projCropland=projectRaster(cropland, crs=aea)
plot(projCropland)
#project a vector
# Datasets
# * CSV table of (fictionalized) brown bear sightings in Alaska, each
# containing an arbitrary ID and spatial location specified as a
# lat-lon coordinate pair.
# * Polygon shapefile containing the boundaries of US National Parks
# greater than 100,000 acres in size.
bears <- read.csv("Data/bear-sightings.csv")
coordinates(bears) <- c("longitude", "latitude")
# read in National Parks polygons
parks <- readOGR("Data", "10m_us_parks_area")
## OGR data source with driver: ESRI Shapefile
## Source: "Data", layer: "10m_us_parks_area"
## with 61 features
## It has 8 fields
# tell R that bear coordinates are in the same lat/lon reference system as the parks data
proj4string(bears) <- proj4string(parks)
# combine is.na() with over() to do the containment test; note that we
# need to "demote" parks to a SpatialPolygons object first
inside.park <- !is.na(over(bears, as(parks, "SpatialPolygons")))
# calculate what fraction of sightings were inside a park
mean(inside.park)
## [1] 0.1720648
## [1] 0.1720648
# determine which park contains each sighting and store the park name as an attribute of the bears data
bears$park <- over(bears, parks)$Unit_Name
# draw a map big enough to encompass all points, then add in park boundaries superimposed upon a map of the United States
plot(bears)
map("world", region="usa", add=TRUE)
plot(parks, border="green", add=TRUE)
legend("topright", cex=0.85, c("Bear in park", "Bear not in park", "Park boundary"), pch=c(16, 1, NA), lty=c(NA, NA, 1), col=c("red", "grey", "green"), bty="n")
title(expression(paste(italic("Ursus arctos"), " sightings with respect to national parks")))
# plot bear points with separate colors inside and outside of parks
points(bears[!inside.park, ], pch=1, col="gray")
points(bears[inside.park, ], pch=16, col="red")
# write the augmented bears dataset to CSV
write.csv(bears, "bears-by-park.csv", row.names=FALSE)
# ...or create a shapefile from the points
writeOGR(bears, ".", "bears-by-park", driver="ESRI Shapefile")
tmin=getData('worldclim', var='tmin', res=10)
# Raster calculator
diff=tmin$tmin1 - tmin$tmin10
## the following code is faster for large datasets.
overlay(tmin$tmin1, tmin$tmin10, fun=function(x,y){return (x-y)})
## class : RasterLayer
## dimensions : 900, 2160, 1944000 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : -343, 102 (min, max)
elevation <- getData("alt", country = "ESP")
slope <- terrain(elevation, opt = "slope")
aspect <- terrain(elevation, opt = "aspect")
hill <- hillShade(slope, aspect, 40, 270)
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Spain")
plot(elevation, col = rainbow(25, alpha = 0.35), add = TRUE)
#contours
contour(elevation)
#crop raster
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Spain")
plot(elevation, col = rainbow(25, alpha = 0.35), add = TRUE)
extent=drawExtent()
cropElev <- crop(elevation, extent)
plot(cropElev)